This script provides a mean of supersaturation time for both species

#load the various libraries
library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(mmtable2)
library(gt)
library(purrr)
library(stringr)
library(tidyr)

Load the data, assign some of the variables as factors

ek_irrad_data <- read.csv("input_data/mean_supersaturation_by_period.csv")

# assign run as a factor
ek_irrad_data$Run <- as.factor(ek_irrad_data$Run)

#assign temperature as a factor
ek_irrad_data$Temperature <- as.factor(ek_irrad_data$Temp...C.)

#assigns treatment as characters from integers then to factors
ek_irrad_data$Treatment <- as.factor(as.character(ek_irrad_data$Treatment))

# make a new column for delta Ek
ek_irrad_data$deltaEk <- (((ek_irrad_data$day9_ek - ek_irrad_data$ek.1)/ ek_irrad_data$ek.1) * 100)

ULVA ____________________________________________________________

Subset the data by species and remove treatment 2.5 because this had problems with tissue sloughing at/near full moon

ulva <- subset(ek_irrad_data, Species == "ul" & Treatment != 2.5)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol" 
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol" 
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol" 
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"

Make a histogram of the data for ulva

hist(ulva$supersat_avg, main = paste("Ulva lactuca"), col = "olivedrab3", labels = TRUE)

ulva %>% ggplot(aes(supersat_avg)) +
        geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
        theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

Run model without interaction between the treatments and temperature - supersat_avg is in minutes Take RLC.Order our of the random effects because causing problems of singularity. R2 is same with or without

hsat_model_ulva <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = ulva)

Check residual plots

hist(resid(hsat_model_ulva))

plot(resid(hsat_model_ulva) ~ fitted(hsat_model_ulva))

qqnorm(resid(hsat_model_ulva))
qqline(resid(hsat_model_ulva))

Check the performance of the model, get R2, summary and print table and plot

performance::check_model(hsat_model_ulva)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(hsat_model_ulva)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m       R2c
## [1,] 0.1963057 0.8597237
summary(hsat_model_ulva)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
##    Data: ulva
## 
## REML criterion at convergence: 2331.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1705 -0.3833  0.0753  0.4654  2.8339 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Plant.ID (Intercept)  609.2   24.68   
##  Run      (Intercept) 2724.1   52.19   
##  Residual              704.8   26.55   
## Number of obs: 240, groups:  Plant.ID, 96; Run, 7
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)    472.758     37.549   5.055  12.591 5.21e-05 ***
## Treatment1     -50.383     44.311   5.002  -1.137    0.307    
## Treatment2     -47.642     44.311   5.002  -1.075    0.331    
## Treatment3     -80.691     44.311   5.002  -1.821    0.128    
## Treatment4     -88.422     44.311   5.002  -1.995    0.103    
## Temperature27    7.294      7.839  75.122   0.930    0.355    
## Temperature30    4.850      7.839  75.122   0.619    0.538    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.835                                   
## Treatment2  -0.835  0.993                            
## Treatment3  -0.835  0.993  0.993                     
## Treatment4  -0.835  0.993  0.993  0.993              
## Temperatr27 -0.104  0.000  0.000  0.000  0.000       
## Temperatr30 -0.104  0.000  0.000  0.000  0.000  0.500
plot(allEffects(hsat_model_ulva))

tab_model(hsat_model_ulva, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  supersat_avg
Predictors Estimates std. Error CI Statistic p df
(Intercept) 472.76 37.55 398.77 – 546.74 12.59 <0.001 230.00
Treatment [1] -50.38 44.31 -137.69 – 36.93 -1.14 0.257 230.00
Treatment [2] -47.64 44.31 -134.95 – 39.67 -1.08 0.283 230.00
Treatment [3] -80.69 44.31 -168.00 – 6.62 -1.82 0.070 230.00
Treatment [4] -88.42 44.31 -175.73 – -1.11 -2.00 0.047 230.00
Temperature [27] 7.29 7.84 -8.15 – 22.74 0.93 0.353 230.00
Temperature [30] 4.85 7.84 -10.60 – 20.30 0.62 0.537 230.00
Random Effects
σ2 704.81
τ00 Plant.ID 609.21
τ00 Run 2724.09
ICC 0.83
N Run 7
N Plant.ID 96
Observations 240
Marginal R2 / Conditional R2 0.196 / 0.860

Construct null model to perform likelihood ratio test REML must be FALSE

ulva_hsat_treatment_null <- lmer(formula = supersat_avg ~ Temperature + (1 | Run) + (1 | Plant.ID) +
                                         (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_hsat_model2 <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +
                                 (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_hsat_treatment_null, ulva_hsat_model2)
## Data: ulva
## Models:
## ulva_hsat_treatment_null: supersat_avg ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_hsat_model2: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                          npar    AIC    BIC  logLik deviance  Chisq Df
## ulva_hsat_treatment_null    7 2461.2 2485.6 -1223.6   2447.2          
## ulva_hsat_model2           11 2396.4 2434.7 -1187.2   2374.4 72.754  4
##                          Pr(>Chisq)    
## ulva_hsat_treatment_null               
## ulva_hsat_model2          5.948e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ulva_hsat_temperature_null <- lmer(formula = supersat_avg ~ Treatment + (1 | Run) + (1 | Plant.ID) +
                                           (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_hsat_model3 <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +
                                 (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_hsat_temperature_null, ulva_hsat_model3)
## Data: ulva
## Models:
## ulva_hsat_temperature_null: supersat_avg ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_hsat_model3: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                            npar    AIC    BIC  logLik deviance  Chisq Df
## ulva_hsat_temperature_null    9 2393.4 2424.7 -1187.7   2375.4          
## ulva_hsat_model3             11 2396.4 2434.7 -1187.2   2374.4 0.9628  2
##                            Pr(>Chisq)
## ulva_hsat_temperature_null           
## ulva_hsat_model3               0.6179

Make plots of the data

ulva %>% ggplot(aes(treatment_graph, supersat_avg)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
        labs(x="Treatment (salinty/nitrate)", y= "Hsat Time (min)", title= "A", subtitle = "Ulva lactuca") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(140, 620) + stat_mean() + 
        geom_hline(yintercept=400, color = "red", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#295102", "#7CB950", "#BDE269")) +
        theme_bw() +
        theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

ulva %>% ggplot(aes(treatment_graph, deltaEk)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = FALSE) + 
        labs(x="Treatment (salinty/nitrate)", y= "delta Ek (%)", title= "A", subtitle = "Ulva lactuca") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(-80, 120) + stat_mean() + 
        geom_hline(yintercept=0, color = "orange", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
        theme_bw() +
        theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))

Summarize the means for variables of interest

ulva %>% group_by(Treatment) %>% summarise_at(vars(supersat_avg), list(mean = mean))
## # A tibble: 5 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          477.
## 2 1          432.
## 3 2          435.
## 4 3          402.
## 5 4          394.
ulva %>% group_by(Run) %>% summarise_at(vars(supersat_avg), list(mean = mean))
## # A tibble: 7 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      461.
## 2 2      465.
## 3 3      422.
## 4 3.5    352.
## 5 4      349.
## 6 6      452.
## 7 6.5    501.
ulva %>% group_by(Run) %>% summarise_at(vars(day_length_avg), list(mean = mean))
## # A tibble: 7 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      654.
## 2 2      643.
## 3 3      626.
## 4 3.5    621.
## 5 4      611.
## 6 6      626.
## 7 6.5    633.
ulva %>% group_by(Treatment) %>% summarise_at(vars(deltaEk), list(mean = mean))
## # A tibble: 5 × 2
##   Treatment   mean
##   <fct>      <dbl>
## 1 0         -27.8 
## 2 1         -22.8 
## 3 2         -22.9 
## 4 3           1.83
## 5 4           4.61

Add growth rate to this dataset. Open growth dataset, prepare data for insertion

growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)
growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100
gr_ulva <- subset(growth_rate, Species == "Ul" & treatment != 2.5)
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)

Plot a regression between the photosynthetic independent variables of interest and growth rate

ulva_growth_hsat_graph <- ggplot(ulva, aes(x=supersat_avg, y=growth_rate)) + 
        geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
        geom_smooth(method = "lm", col = "black") + theme_bw() + 
        labs(title = "Ulva lactuca Hsat vs Growth Rate", x = "Hsat time (minutes)", 
             y = "growth rate (%)") + stat_regline_equation(label.x = 240, label.y = 165) + stat_cor()
ulva_growth_hsat_graph
## `geom_smooth()` using formula = 'y ~ x'

HYPNEA_____________________________________________________________________________________________________________

hypnea <- subset(ek_irrad_data, Species == "hm" & day1_rlc_time != "11:34:05")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol" 
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol" 
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol" 
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"

Make a histogram for Hypnea

hist(hypnea$supersat_avg, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)

hypnea %>% ggplot(aes(supersat_avg)) +
        geom_histogram(binwidth=5, fill = "#a8325e", color = "black", size = 0.25, alpha = 0.85) +
        theme_bw()

Run model for hsat hypnea

hsat_model_hypnea <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = hypnea)

Make residual plots of the data

hist(resid(hsat_model_hypnea))

plot(resid(hsat_model_hypnea) ~ fitted(hsat_model_hypnea))

qqnorm(resid(hsat_model_hypnea))
qqline(resid(hsat_model_hypnea))

Check the performance of the model, get r2, plot model and print table of the data

performance::check_model(hsat_model_hypnea)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(hsat_model_hypnea)
##           R2m       R2c
## [1,] 0.144871 0.7606678
summary(hsat_model_hypnea)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
##    Data: hypnea
## 
## REML criterion at convergence: 2894.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3056 -0.4446  0.0506  0.5483  2.1875 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Plant.ID (Intercept)  790.2   28.11   
##  Run      (Intercept) 2049.5   45.27   
##  Residual             1103.6   33.22   
## Number of obs: 286, groups:  Plant.ID, 144; Run, 8
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)    396.461     32.941   5.042  12.035 6.62e-05 ***
## Treatment1     -49.575     38.936   5.020  -1.273 0.258708    
## Treatment2     -35.802     38.936   5.020  -0.920 0.399862    
## Treatment2.5   -37.009     56.153   4.730  -0.659 0.540567    
## Treatment3     -51.422     38.936   5.020  -1.321 0.243612    
## Treatment4     -69.064     38.951   5.027  -1.773 0.136083    
## Temperature27   29.379      7.917 131.257   3.711 0.000304 ***
## Temperature30   32.503      7.927 132.169   4.100 7.17e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.830                                          
## Treatment2  -0.830  0.985                                   
## Treatmnt2.5 -0.575  0.487  0.487                            
## Treatment3  -0.830  0.985  0.985  0.487                     
## Treatment4  -0.830  0.984  0.984  0.487  0.984              
## Temperatr27 -0.120  0.000  0.000  0.000  0.000  0.000       
## Temperatr30 -0.120  0.000  0.000  0.000  0.000  0.001  0.499
tab_model(hsat_model_hypnea, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  supersat_avg
Predictors Estimates std. Error CI Statistic p df
(Intercept) 396.46 32.94 331.61 – 461.31 12.04 <0.001 275.00
Treatment [1] -49.58 38.94 -126.23 – 27.08 -1.27 0.204 275.00
Treatment [2] -35.80 38.94 -112.45 – 40.85 -0.92 0.359 275.00
Treatment [2.5] -37.01 56.15 -147.55 – 73.53 -0.66 0.510 275.00
Treatment [3] -51.42 38.94 -128.07 – 25.23 -1.32 0.188 275.00
Treatment [4] -69.06 38.95 -145.74 – 7.62 -1.77 0.077 275.00
Temperature [27] 29.38 7.92 13.79 – 44.97 3.71 <0.001 275.00
Temperature [30] 32.50 7.93 16.90 – 48.11 4.10 <0.001 275.00
Random Effects
σ2 1103.65
τ00 Plant.ID 790.18
τ00 Run 2049.48
ICC 0.72
N Run 8
N Plant.ID 144
Observations 286
Marginal R2 / Conditional R2 0.145 / 0.761 
plot(allEffects(hsat_model_hypnea))

Construct null model to perform likelihood ratio test REML must be FALSE

hyp_hsat_treatment_null <- lmer(formula = supersat_avg ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hyp_hsat_model2 <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hyp_hsat_treatment_null, hyp_hsat_model2)
## Data: hypnea
## Models:
## hyp_hsat_treatment_null: supersat_avg ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hyp_hsat_model2: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                         npar    AIC    BIC  logLik deviance  Chisq Df
## hyp_hsat_treatment_null    7 2987.0 3012.6 -1486.5   2973.0          
## hyp_hsat_model2           12 2972.1 3015.9 -1474.0   2948.1 24.914  5
##                         Pr(>Chisq)    
## hyp_hsat_treatment_null               
## hyp_hsat_model2          0.0001448 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hyp_hsat_temperature_null <- lmer(formula = supersat_avg ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hyp_hsat_model3 <- lmer(formula = supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hyp_hsat_temperature_null, hyp_hsat_model3)
## Data: hypnea
## Models:
## hyp_hsat_temperature_null: supersat_avg ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hyp_hsat_model3: supersat_avg ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                           npar    AIC    BIC  logLik deviance  Chisq Df
## hyp_hsat_temperature_null   10 2985.0 3021.5 -1482.5   2965.0          
## hyp_hsat_model3             12 2972.1 3015.9 -1474.0   2948.1 16.894  2
##                           Pr(>Chisq)    
## hyp_hsat_temperature_null               
## hyp_hsat_model3            0.0002145 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Make plots

hypnea %>% ggplot(aes(treatment_graph, supersat_avg)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
        labs(x="Treatment (salinty/nitrate)", y= "Hsat Time (min)", title= "B", subtitle = "Hypnea musciformis") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(140, 620) + stat_mean() + 
        geom_hline(yintercept=400, color = "red", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#9C0627", "#BB589F", "#F4B4E2")) +
        theme_bw() +
        theme(legend.position = c(0.88,0.88), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

Plot deltaEk for Hypnea

hypnea %>% ggplot(aes(treatment_graph, deltaEk)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = FALSE) + 
        labs(x="Treatment (salinty/nitrate)", y= "delta Ek (%)", title= "B", subtitle = "Hypnea musciformis") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(-100, 250) + stat_mean() + 
        geom_hline(yintercept=0, color = "orange", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
        theme_bw() +
        theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))

Prep data for linear regression with growth. For Hypnea remove hm6-4 on 11/12 that had no d9 RLC (final weight 0.1017) and hm6-4 on 10/09/21 because it was white and also looked dead

gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)

Plot a regression between the photosynthetic independent variables of interest and growth rate

hypnea_growth_irrad_ek_graph <- ggplot(hypnea, aes(x=supersat_avg, y=growth_rate)) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
        geom_smooth(method = "lm", col = "black") + theme_bw() + 
        labs(title = "B", subtitle = "Hypnea musciformis", x = "Hsat Time (min)", y = "growth rate (%)") + 
        stat_regline_equation(label.x = 400, label.y = 150) + stat_cor(label.x = 400, label.y = 140) +
        theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
hypnea_growth_irrad_ek_graph
## `geom_smooth()` using formula = 'y ~ x'

Summarize the means for rETRmax

hypnea %>% group_by(Treatment) %>% summarise_at(vars(supersat_avg), list(mean = mean))
## # A tibble: 6 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          417.
## 2 1          374.
## 3 2          388.
## 4 2.5        380.
## 5 3          373.
## 6 4          354.
hypnea %>% group_by(Run) %>% summarise_at(vars(supersat_avg), list(mean = mean))
## # A tibble: 8 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      407.
## 2 2      427.
## 3 3      360.
## 4 3.5    312.
## 5 4      320.
## 6 7      380.
## 7 8      406.
## 8 9      428.
hypnea %>% group_by(Run) %>% summarise_at(vars(day_length_avg), list(mean = mean))
## # A tibble: 8 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      654.
## 2 2      645.
## 3 3      626.
## 4 3.5    621.
## 5 4      611.
## 6 7      688.
## 7 8      635.
## 8 9      636.
hypnea %>% group_by(Treatment) %>% summarise_at(vars(deltaEk), list(mean = mean))
## # A tibble: 6 × 2
##   Treatment   mean
##   <fct>      <dbl>
## 1 0           9.46
## 2 1         -37.4 
## 3 2         -33.2 
## 4 2.5         7.35
## 5 3         -17.1 
## 6 4         -21.1